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Abstract 

We report next-to-leading order perturbative QCD predictions of 4 jet event shape 
variables for the process e"'"e~ — > 4 jets obtained using the general purpose Monte Carlo 
EERAD2. This program is based on the known 'squared' one loop matrix elements for 
the virtual 7* — > 4 parton contribution and squared matrix elements for 5 parton 
production. To combine the two distinct final states numerically we present a hybrid 
of the commonly used subtraction and slicing schemes based on the colour antenna 
structure of the final state which can be readily applied to other processes. We have 
checked that the numerical results obtained with EERAD2 are consistent with next- 
to- leading order estimates of the distributions of previously determined four jet-like 
event variables. We also report the next-to-leading order scale independent coefficients 
for some previously uncalculated observables; the light hemisphere mass, narrow jet 
broadening, Aplanarity and the 4 jet transition variables with respect to the JADE and 
Geneva jet finding algorithms. For each of these observables, the next-to-leading order 
corrections calculated at the physical scale significantly increase the rate compared to 
leading order (the K factor is approximately 1.5 - 2). With the exception of the 4 jet 
transition variables, the published DELPHI data lies well above the ) predictions. 
The renormalisation scale uncertainty is still large and in most cases the data prefers 
a scale significantly smaller than the physical scale. This situation is reminiscent of 
that for three jet shape variables, and should be improved by the inclusion of power 
corrections and resummation of large infrared logarithms. 



1 Introduction 



Electron-positron colliders, in particular those at both LEP and SLAC, have provided much 
precision data with which to probe the structure of QCD. This is particularly valuable data 
because the strong interactions occur only in the final state and are not entangled with 
the parton density functions associated with beams of hadrons. In addition to measuring 
multi-jet production rates, more specific information about the topology of the events can be 
extracted. To this end, many variables have been introduced which characterize the hadronic 
structure of an event. For example, we can ask how planar or how collimated an event is. 
In general, a variable is described as n jet-like if it vanishes for a final state configuration of 
n — 1 hadrons. With the precision data from LEP and SLC, experimental distributions for 
such event shape variables have been studied and have been compared where possible with 
theoretical calculations. 

Generally speaking, leading order (LO) predictions successfully predict the general fea- 
tures of distributions, but can be improved by resumming kinematically-dominant loga- 
rithms, by including more perturbative information or both. A next-to-leading order (NLO) 
treatment of three-jet like variables was first performed in [||, Q and systematically com- 
pleted in . Armed with such calculations, one can extract a value for the strong coupling 
tts either directly from the event shape distributions 0] or from the energy dependence of 
their average value [Q. Alternatively, one can study the group parameters of the gauge the- 
ory of the strong interactions 0, though for three jet observables, the gluon-gluon coupling 
(proportional to Ca) occurs first at NLO. 

Recent calculations of the relevant one-loop four parton matrix elements for 7* — 4 par- 
tons 1^ and e^e~ — >■ 4 partons [Q, together with the known tree-level five parton matrix 
elements have enabled the phenomenology of four-jet production to be be studied at 
next-to-leading order. So far, two groups have used these matrix elements to construct 
general purpose Monte Carlo programs for four jet-like quantities MENLO PARC [|T0| and 
DEBRECEN [|ll|, which have been used to measure the four jet fraction i?4 and a variety 



of event shape distributions [|TT]. We note that four jet production is sensitive to the casimir 
structure of QCD || and four jet events may be used to constrain the allowed values 
of Cp, Ca and Tr |]T3[ by examining the angles between the jets One may also place 



next-to-leading order bounds on the possible presence of the elusive light gluino present in 
many supersymmetric models 



Any such Monte Carlo program must suitably combine the virtual four parton and real 
five parton pieces of the cross section. That is, 

c^NLo = j da, + j da, = j dPS^ |A^4|' + / dPS, \M,\\ (1.1) 

where the n-parton contributions dcr^ are integrals over the n-parton phase-space dPSn- 
Both integrals are separately divergent, and contain both infrared and ultraviolet singular- 
ities. The ultraviolet poles are removed by renormalisation, however the soft and collinear 
infrared poles are only cancelled when the virtual graphs are combined with the brem- 
strahlung process and one must provide a means to cancel the infrared singularities caused 
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by two particles becoming collinear or one soft. Although the cancellation of infrared poles 
can be done analytically for simple processes, for complicated processes, it is necessary to 
resort to numerical techniques. A variety of methods - known as subtraction 16, l7 



slicing p|, |T8[ and hybrid subtraction [|T9l - are in general use, and variations of the sub- 



traction formalism have been implemented in |T0, |Tl|. Here we report on results obtained 



using a third numerical implementation of these matrix elements to compute generic infrared 



safe four jet observables. We call this program EERAD2 and it is based on the 'squared' 
one-loop matrix elements for 7* — > 4 partons of together with squared tree level matrix 
elements for 7* — > 5 partons. To isolate and cancel the infrared singularities, we use a hybrid 
approach that contains elements of both slicing and subtraction methods. In particular, we 
use an antenna factorisation where two (colour connected) hard partons radiate a third that 
may be unresolved. This is necessarily a rather technical subject and it is detailed in the 
appendix. 

For the bulk of the paper we are more concerned with the phenomenology of four jet-like 
shape variables and, in particular, we extend the set of variables that have been computed 
at 0(0:3). [| To be precise, we present next-to-leading order coefficents for the differential 
distributions of the narrow jet broadening, light hemisphere mass, Aplanarity and the jet 
transition variable for the JADE and Geneva jet algorithms. To compare with previous 
results, we also consider the thrust minor, D parameter and the jet transition variable for 
the Durham jet algorithm as well as the four jet rate as a function of the jet resolution 
parameter ?/cut- 

In section 2, we give the definitions of the relevant four jet shape variables and review 
the structure of the perturbative predictions. Section 3 shows the consistency of our pro- 
gram, by reproducing the known four-jet fraction and D parameter distributions. The main 
results are reported in section 4 and comparisons with experimental data follow in section 5, 
where we also discuss how we might optimize our perturbative input by choice of a suitable 
renormalization scale. Conclusions are summarized in section 6. Finally, a more detailed 
discussion of our method for cancellation of the infrared singularities is reserved for the 
Appendix. 



Some of the results presented here have been reported in [ElJ 
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2 Four jet event shapes 



The sorts of variables we are interested in are four jet-like, since they can only be non-zero 
for final states in which there are four or more particles. They usually rely on the hadronic 
final state having some volume and, when the event is coplanar, some observables like the 
D parameter are identically zero. 



2.1 Definition of Variables 

In the following definitions, the sums run over all final state particles, k = 1, . . . ,N . Pk is 
the three-momentum of particle k in the cm. frame, with components p^, i = 1,2,3. 



(a) C and D parameters 



We first construct the linear momentum tensor, 

\Pk\ 



e^j = IP; I , (2.1) 



with eigenvalues Aj for i = 1,2, 3. The normalisation is such that X^i = 1- For planar 
events one of the eigenvalues is zero. The C and D parameters are defined by, 

D = 27A1A2A3, (2.2) 

and, 

C = 3(AiA2 + A2A3 + A3Ai). (2.3) 

D can only be non-zero for non-planar four (or more) parton events, while three parton 
events may produce < C < 0.75. Only the region C > 0.75 should be considered 
four jet-like. 



(b) Thrust minor, Tminor ||23 



We first define the thrust, major and minor axes {ni,n2, n^) by, 

T. = max5|ML^, (2.4) 

Efc \Pk\ 

where n2 is constrained by ni ■ ^2 = 0. and ns = ni x n2. 
(c) Light hemisphere mass, M\l s. 

The event is separated into two hemispheres H2 divided by the plane normal to 
the thrust axis ni, as defined above. Particles that satisfy Pi.ni > are assigned to 
hemisphere Hi, while all other particles are in H2. Then, 

E Pk] ■ (2.5) 



Ml 
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Note that this is the common modification of the original definition suggested by 
Clavelh l2l. 



(d) Narrow jet broadening, -Bmin p5[] . 

Using the same division into hemispheres as above, we define. 



i=l,2 





Pk X n\ 




\Pk 





-Dmin — min . {^■0) 



(e) Aplanarity, A ^6 . 



Here we consider the eigenvalues, Ai > A2 > A3 of the quadratic momentum tensor. 



^ij ^ g^M. (2.7) 



Aplanarity is defined by, 

A = hs, (2.8) 
which is clearly zero for planar three particle events. 

(f) Jet transition variable yf. 

The j/f variable denotes the value of the jet resolution parameter ?/cut at which an event 
changes from a four jet event to a three jet event where the jets are defined according 
to algorithm S. We consider three algorithms, the JADE algorithm {S = J) ||27[| , 
the Durham algorithm {S = D) p8| and the Geneva algorithm {S = G) p9|. The 



jet-finding measures for each of these three algorithms are as follows, 

J _ 2EiEj{l -cos%) 

^ _ 2 mm{Ef, E]){1 - cos%) 
yij - I ' 

G _ 8 EjEjjl- cos Oij) 
^'^~9 {E, + E,y ' ^^-^^ 

where the factor of 8/9 in the Geneva algorithm is simply to ensure that the maximum 
value of i/cut that reconstructs three jets from three partons is 1/3 as it is for the other 
two algorithms. When particles combine, there is some ambiguity as to how to add 
the energies and momenta. In all three schemes, we use the E scheme i.e. we merely 
add four momenta, 

pf, =pr+p^ (2.10) 

Other choices such as the EO or P schemes where the cluster is made massless by 
rescaling the momentum or energy give similar results. 



Of these variables, the D, C, Tminor and yf distributions have been studied in 
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2.2 Structure of Perturbative Prediction 



The differential cross-section at centre-of-mass energy ^/s for one of tliese four-jet variables 
(O4) at next-to-leading order is described by two coefficients, Bq^ and Cq^ which represent 
the leading and next-to-leading order perturbative contributions, 

y.--^^r (^)' - {"^J (t) -^o. . c.) , (.„, 

Both Bq^ and Co^ are scale independent and do not depend on the beam energy. However, 
the running coupling is calculated at renormalization scale /i which is commonly chosen 
to be the physical scale, fi = y/s. Compared to the leading order prediction, which decays 
monotonically with increasing /x, the next-to- leading order term reduces the scale dependence 
somewhat through the first coefficient of the beta-function, (3o = {33—2Nf) / 6. For five active 
quark flavours, Po = 3.833. 

The individual coefficients B04 and C04 depend on the number of colours and flavours, 
or equivalently, the group casimirs of the standard model. As such, four jet event shapes 
may be used to simultaneously constrain the strong interaction gauge group as well as the 
strong coupling constant. 



2.3 Scale choice and theoretical uncertainty 

As mentioned above, for hadronic observables in electron-positron annihilation it is common 
to choose the renormalisation scale to be the physical scale /z = ^/s. This choice is motivated 
by naturalness arguments and the fact that choosing a scale far from ^/s introduces large 
logarithms of the form \og{^/y/s) in eq. ( |2.11| ). Since the renormalisation scale is only a 



theoretical construct, it has become common practice to estimate the uncertainty engendered 
by truncating the perturbative expansion by varying /i by factors about the physical scale. 

Other approaches have been considered with rather different scale choices. One such 
approach is to stipulate that the next-to-leading order coefficient vanishes. That is to say 
that the LO and NLO predictions coincide.| This occurs at the scale, 

"^^^ - ^ (-silk) • <-^) 

and is called the scale of fastest apparent convergence (FAC) As we will see, C04 is 

typically 50-100 times larger than B04 so that fj,^^'" may be as small as 10~^-10~^ of the 
physical scale. 

Another common choice is to select the scale where the NLO prediction is most insensitive 



to the choice of scale |3l|. However, this scale is very close to the FAC scale (15% smaller) so 



^In other words, the ratio of NLO to LO predictions commonly called the K-factor is identically unity. 
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that the Principle of Minimal Sensitivity scale (PMS) predictions he very close to the FAC 
scale prediction. 

While even higher order corrections remain uncalculated, varying the renormalisation 
scale can only give a crude indication of the theoretical uncertainty. Therefore, in an attempt 
to make a fair estimate of the theoretical uncertainty on the NLO prediction we will show 
both the physical scale and FAC scale predictions. 

2.4 Infrared behaviour 

Four jet event shapes typically depend on the event having some volume and not lying entirely 
in a plane. Typical hadronic events contain more than 20 hadrons and it is extremely unlikely 
that the value of any event shape is precisely zero for any experimental event. However, in 
a LO or NLO fixed order parton calculation, there only four or five partons present in the 
final state and, when one or more are soft, the calculated O4 may approach zero. In such 
circumstances, soft gluon singularities cause the fixed order prediction to become wildly 
unstable and grow logarithmically. In the small O4 limits, the perturbative coefficients have 
the following form, 

B04 A32L^ + A22L^ + AuL + A02, 

Co, ^ A^^L'' + A^^L^ + A^^L"" + A2^L^ + Ai^L + Aoz, (2.13) 

where L = log(l/04) and Anm are (as yet) undetermined coefficients. Whenever L is suffi- 
ciently large, resummation effects will be important .0 In comparing with data, we therefore 
choose to make a cut on the size of O4 which is typically in the range 0.001 — 0.01, since for 
such small values of O4 we do not trust the NLO prediction. In comparing with the DELPHI 
data, this cut will usually be the lower edge of the second data bin. 



■^Whether the coefficients exponentiate and can be resummcd will depend on the observable. 
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3 Comparison with existing results 



3.1 Four jet rates 

As a check of the numerical results, Table |I| shows the predictions for each of the three 
Monte Carlo programs for the four jet rate for three jet clustering algorithms; the Jade-EO, 
Durham-E and Geneva-E |29] algorithms. We show results with as{Mz) = 0.118 for 



three values of the jet resolution parameter i/cut- There is good agreement with the results 
from the other two calculations. 



Algorithm 


Vcut 


MENLO PARC 


DEBRECEN 


EERAD2 


Durham 


0.005 
0.01 
0.03 


(1.04 ±0.02) ■ 10"^ 
(4.70 ±0.06) ■ 10-2 
(6.82 ±0.08) ■ 10-3 


(1.05 ±0.01) ■ 10-1 
(4.66 ±0.02) ■ 10-2 
(6.87 ±0.04) • 10-3 


(1.05 ±0.01) ■ 10-1 
(4.65 ±0.02) ■ 10-2 
(6.86 ±0.03) • 10-3 


Geneva 


0.02 
0.03 
0.05 


(2.56 ±0.06) ■ 10-1 
(1.71 ±0.03) ■ 10-1 
(8.58 ±0.15) ■ 10-2 


(2.63 ±0.06) ■ 10-1 
(1.75 ±0.03) ■ 10-1 
(8.37 ±0.12) ■ 10-2 


(2.61 ±0.05) ■ 10-1 
(1.72 ±0.03) ■ 10-1 
(8.50 ±0.06) ■ 10-2 


JADE-EO 


0.005 
0.01 
0.03 


(3.79 ±0.08) • 10-1 
(1.88 ±0.03) • 10-1 
(3.46 ±0.05) • 10-2 


(3.88 ±0.07) ■ 10-1 
(1.92 ±0.01) • 10-1 
(3.37 ±0.01) ■ 10-2 


(3.87 ±0.03) ■ 10-1 
(1.93 ±0.01) ■ 10-1 
(3.35 ±0.01) ■ 10-2 



Table 1: The four-jet fraction as calculated by MENLO PARC, DEBRECEN and EERAD2, for the 
different jet recombination schemes and varying ?/cut- The rate is normalized by the 0{as) 
total hadronic cross-section, cThad = c"o (1 ± ^s/tt). 



3.2 Shape variables 



As mentioned earlier, Nagy and Trocsanyi [O have computed Cd with their Monte Carlo 



DEBRECEN. In Table ^ we show the leading and next-to- leading order coefficients and Cd 
calculated by EERAD2, together with the DEBRECEN result. The two calculations are clearly 
consistent with one another, with the quoted errors overlapping in almost all cases. The 
errors from EERAD2 are of the order of 2% in each bin, except in the tail of the distribution 
where the errors rise as high as 10%. The infrared enhancement of the distribution described 
in section |2.4| means that the Monte Carlo procedure favours the phase space region corre- 



sponding to small values of the D parameter, so that the large D tail suffers larger errors. 
In fact Cd drops by four orders of magnitude over the kinematic range of the observable 
so it is necessary to use importance sampling with respect to the observable distribution to 
ensure sufficient Monte Carlo points are produced in the high D region. This is also true for 
all of the other shape variables. 

In addition, Nagy and Trocsanyi have also presented results for the next-to-leading order 
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coefiicents for thrust minor Tminor and the jet transition variable in the Durham scheme 
EH- Although we do not present a detailed comparison here, we note that the agreement 
is qualitatively the same as discussed for the D parameter above. We find that the distri- 
butions extend beyond the range of coefficents presented in with non-zero coefficients 
for bins in the ranges 0.5 < Tminor < 0.58 and 0.125 < < 0.17. 



D 


Bd 


c 







DEBRECEN 


0200 


('Q 7Q -1- n ni 
yo. 1 y ic u.ui ) 


• iU 


fl 47 ± 


(]()) 


• 10^ 

X w 


(-] ns -1- n r\R\ 

l^i.Uo It U.UO ) ■ 


iU 


0600 


yz.oz zsz u.ui J 


■ iU 


(1.25 ± 


01) 


■ 10^ 


l^i.Z^t It U.UZ J • 


iU 


1000 


yL.^o It u.ui ) 


• iU 


f8 69 ± 


04) 


• 10^ 

X w 


yo.Ot) It U. i/J ■ 


iU 


1400 


/ -\ n/i 4- n ni 

(^i.U4 It U.Ui ) 


■ iU 


(6 39 ± 


03^ 


• 10^ 

X w 


((\ OA 4- r\ T 0\ 
\K).Z^ It U. iZJ ■ 


iU 


1800 


y 1 .oo It U.U4I: J 


ini 

• iU 


f4 89 ± 


03^ 


• 10^ 

X w 


y^.ifi) It u. i i J ■ 


iU 


2200 


\o.O 1 It U.UO ) 


• iU 


(3.88 ± 


03) 


• 10^ 


yo.oo It U.UO J ■ 


iU 


2600 


^^41.00 It U.U ( J 


■ iU 


(3.04 ± 


03) 


• 10^ 


yz.ilo It U.UO J ■ 


iU 


3000 

\J • KJ\J\J\J 


yo. ( It u.u ( J 


1 ni 

• iU 


(2.51 ± 


04) 


• 10^ 


(0 KO 4- C\ C\K\ 
yZ.oZ It U.UO J • 


iU 


3400 


yo.U 1 It U.UO ) 


■ iU 


(2.02 ± 


03") 

WO 1 


• 10^ 

X w 


l^i.y^i: It U.UO ) • 


iU 


3800 


I Z.^l It U.UO I 


• iU 


(1.61 ± 


03) 


• 10^ 


^^i.oy It u.u^ j • 


1 n3 

iU 


0.4200 


{^ Q7 -U n HA') 
^^i.y ( It u.u^ ) 


■ ±u 


(1.37 ± 


02) 


■ 10^ 


yL.O i It U.UO J ■ 


1 n3 

iU 


0.4600 


(1 56 ± 03) 


■ 10^ 


(1.09 ± 


01) 


• 10^ 


fl 06 ± 03) ■ 

I X . W W — 1 — \J •\JkJ j 


10^ 


0.5000 


(1.32 ±0.01) 


■10^ 


(8.97 ±0 


14) 


■ 102 


(8.72 ±0.19) ■ 


102 


0.5400 


(1.05 ±0.02) 


•101 


(7.12 ±0 


15) 


■ 102 


(7.11 ±0.16) ■ 


102 


0.5800 


(8.46 ±0.16) 


■10° 


(5.79 ±0 


12) 


• 102 


(5.68 ±0.14) ■ 


102 


0.6200 


(6.60 ±0.16) 


■10° 


(4.55 ± 


09) 


■ 102 


(4.46 ±0.21) ■ 


102 


0.6600 


(5.32 ±0.13) 


•10° 


(3.58 ±0 


07) 


• 102 


(3.52 ±0.11) ■ 


102 


0.7000 


(3.99 ±0.09) 


•10° 


(2.80 ±0 


09) 


• 102 


(2.74 ±0.09) ■ 


102 


0.7400 


(3.06 ±0.05) 


■10° 


(2.05 ±0 


08) 


■102 


(2.08 ±0.08) ■ 


102 


0.7800 


(2.26 ±0.04) 


■ 10° 


(1.58 ±0 


04) 


■ 102 


(1.54 ±0.06) ■ 


102 


0.8200 


(1.54 ±0.04) 


■10° 


(1.05 ±0 


03) 


■102 


(1.03 ±0.04) ■ 


102 


0.8600 


(9.72 ±0.21) 


10-1 


(6.72 ± 


29) 


■101 


(6.66 ±0.31) ■ 


101 


0.9000 


(5.63 ±0.16) 


10-1 


(3.85 ±0 


17) 


■101 


(3.89 ±0.20) ■ 


101 


0.9400 


(2.62 ±0.07) 


10-1 


(1.71 ±0 


10) 


■lOi 


(1.71 ±0.19) ■ 


lOi 


0.9800 


(5.34 ±0.11) 


10-2 


(3.15 ±0 


27) 


■ 10° 


(2.60 ± 1.30) ■ 


10° 



Table 2: The leading and next-to-leading order coefficients for the D parameter. The NLO 



coefficient predicted by Nagy and Trocsanyi Monte Carlo DEBRECEN |TT| is also shown. 
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4 New results 



In this section we extend the analysis of 4 jet-like event shape observables already found 
in the literature by reporting the leading and next-to-leading order coefficients for the light 
hemisphere mass, the narrow hemisphere broadening, Aplanarity and the jet transition vari- 
able in both the JADE and Geneva schemes, and yf. In particular, we examine the 
relative sizes of the two terms by inspecting the K factor (at the physical scale) for each 
variable across the allowed kinematic range of the distributions. 

For all the variables presented in this section, we must be careful to differentiate between 
the true behaviour of the distribution as the observable tends to zero and the behaviour in 
fixed order perturbation theory. Each of the observables should have a smooth behaviour 
as O4 ^ rather than the divergent behaviour exhibited by the coefficients according to 
equation |2.13 . To recover a smooth result in this limit it is necessary to resum powers of 



log(l/04) where possible, a procedure which has been performed already for many 3 jet-like 
variables [0^, 



4.1 Light Hemisphere Mass 



M'l/s 


BMf/s 




0.0150 


(3.23 ±0.08) 


■10^ 


(1.41 ±0.01) 


■ 10^ 


0.0250 


(1.88 ±0.02) 


■102 


(8.85 ±0.10) 


■ 103 


0.0350 


(1.25 ±0.02) 


■ 10^ 


(5.97±0.11) 


■ 103 


0.0450 


(8.52 ±0.10) 


■W 


(4.14 ±0.08) 


■ 103 


0.0550 


(5.97 ±0.06) 




(3.04 ±0.04) 


■ 103 


0.0650 


(4.20 ±0.09) 




(2.15 ±0.05) 


■ 103 


0.0750 


(3.02 ±0.07) 




(1.58 ±0.05) 


■ 103 


0.0850 


(2.13 ±0.03) 


■101 


(1.11 ±0.02) 


■ 103 


0.0950 


(1.39 ±0.04) 


■10^ 


(7.66 ±0.23) 


■102 


0.1050 


(8.75 ±0.20) 


■10° 


(4.97 ±0.17) 


■102 


0.1150 


(5.18 ±0.13) 


■10° 


(3.27 ±0.07) 


■102 


0.1250 


(2.59 ±0.12) 


■10° 


(1.66 ±0.07) 


■102 


0.1350 


(8.97 ±0.35) ■ 


10-1 


(6.61 ±0.41) 


■101 


0.1450 


(2.49 ±0.13) ■ 


10-1 


(1.79 ±0.09) 


■101 


0.1550 


(5.00 ±0.27) ■ 


10-2 


(3.75 ±0.26) 


■10° 


0.1650 


(1.46 ±0.21) ■ 


10-3 


(2.30 ±0.37) 


10-1 



Table 3: The leading and next-to-leading order coefficients for the light jet mass M'l/s. 



As defined before, the light hemisphere mass is the smaller invariant mass of the two 
hemispheres formed by separating the event by a plane normal to the thrust axis. The NLO 
coefficient C^i/s evaluted at the physical scale fi = ^/s together with the LO term is given 
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in Table |^. The errors are estimates from the numerical program and are typically 2-3% for 
each entry. As with the previously known results on four jet event shapes, the NLO terms 
are significantly larger than the LO term. Here, we see that 0^2/^ is typically 50 times 
larger than Bjy.f2/g so that even when the additional factor of as/27r is restored, the NLO 
correction is large. This is illustrated in Fig. where the K factor defined by, 

is shown for O4 = M|/ s. We see that the K factor increases with the value of the observable, 
rising from 1.8 at small M|/s up to 2.4. This behaviour is similar to that observed for other 
four jet event shapes |[TT|. 




Figure 1: The K-factors defined according to eq. (f4.1|) for four jet event shapes, the light 
hemisphere mass (solid), narrow jet broadening (long-dashed), Aplanarity (dot-dashed) and 
jet transition variables in the JADE (short-dashed) and Geneva (dotted) schemes. Each 
variable has a different kinematic range. 



4.2 Narrow Hemisphere Broadening 

Narrow hemisphere broadening, -Bmin, is defined in a similiar manner to the light hemisphere 
mass. The event is again divided into two hemispheres by the plane normal to the thrust 
axis, but now the momenta transverse to the thrust axis is summed (normalised by the sum 
of absolute momenta) in each hemisphere. The narrow hemisphere is that with the least 
transverse momentum with respect to the thrust axis. Numerical results for this variable 
as calculated by EERAD2 can be found in Table |. As with the light hemisphere mass, the 
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NLO contribution is significant yielding a K factor of roughly 1.7 over most of the kinematic 
range of the variable (see Fig. |^). 
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Table 4: The leading and next-to-leading order coefficients for the narrow jet broadening 



4.3 Aplanarity 

As described earlier, Aplanarity is essentially the smallest eigenvalue of the quadratic mo- 
mentum tensor. The NLO coefficient Ca evaluted at the physical scale /i = y/s together 
with the LO term is given in Table ^. Once again, the NLO terms are significantly larger 
than the LO term Ba and, as can be seen in Fig. |I| gives rise to K factor of roughly 1.5. 



4.4 Jet transition variables 

As previously stated the jet transition variable yf describes the scale where two jets merge, 
thereby changing a four jet event into a three jet event. This is essentially the same as the 
derivative of the four jet rate with respect to the jet resolution parameter ijcnt- However, 
the number of jets in an event is dependent on the jet finding algorithm used to define 



the 'closeness' of particles which is compared with ?/cut- h^ the transition rate for the 
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Table 5: The leading and next-to- leading order coefficients for Aplanarity. 



the Durham jet finding algorithm 28 is given and we have checked that our results are 



consistent with these predictions. Here, we provide results for two other jet algorithms, the 
JADE and Geneva schemes for which the jet finding measures are given in eq. ( |2.9| ). 
We note that the Geneva algorithm enjoys the same benefits as the Durham algorithm in 
that it is also supposed to exponentiate, enabling infrared logarithms to be safely resummed. 
It also ensures that softly radiated gluons are clustered with hard partons unless the angle 
of separation between two soft gluons is much smaller than the angular separation between 
them and a hard parton. 

Our results for the two schemes are given in Tables |^ and |^. As can be seen from the 
tables the NLO coefficients are large, which is again reflected in the large corrections shown 
in Fig. |I|. The K factor for the JADE scheme is roughly 1.8-1.9, but is slightly smaller for 
the Geneva algorithm, typically in the region 1.4-1.6. 
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Table 6: The leading and next-to-leading order coefficients for the jet transition variable in 
the Geneva-E algorithm y^. 
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Table 7: The leading and next-to-leading order coefficients for the jet transition variable in 
the JADE-EO algorithm y(. 
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5 Comparison with experimental data 



Four jet event shape observables have been studied extensively by the four LEP experiments. 
However, the most complete analysis of event shape variables has been carried out by the 
DELPHI collaboration |33|. Here, they study all of the event shape variables discussed in 



section |^. Distributions based on charged particles alone as well as charged and neutral 
particles are presented. In this section, we wish to examine whether or not these event 
shapes can be described by fixed order perturbation theory. As discussed earlier, to avoid 
numerical instabilities in the infrared region where fixed order perturbation theory is no 
longer valid we impose a cut on the smallness of the variable that is generally equal to the 
lower edge of the second bin. More precisely, that is. 
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0.008, 


T ■ 

^ minor 


> 


0.02, 


Ml/s 
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0.01, 




> 


0.01, 


A 


> 


0.005, 




> 


0.002, 


yi 


> 


0.002. 



(5.1) 



The experimental distributions are normalised to the hadronic cross section (rather than 
the Born cross section) and are also not weighted by the observable, but are rather. 



da 




(5.2) 



Throughout, we choose as{Mz) = 0.118 which is consistent with the current world average 



3^ . In each case, the theoretical predictions have been evaluated using bins of the same size 
as in the experiment and therefore appear as histograms in the plots. The data is corrected 
for detector effects, but not for hadronisation effects. 

Figs. 1^, ^ and | show the comparison between the leading order and next-to-leading 
order predictions evaluated at the physical scale /x = y/s = Mz for narrow jet broadening, 
light hemisphere mass and Aplanarity with the published DELPHI data ||3^. We see that 
in all three cases, the LO prediction undershoots the data by a significant factor (about a 
factor of four), and that including the NLO correction improves the situation but still gives 
a rate that is much lower than the data. However, the NLO prediction still contains a large 
renormalisation scale uncertainty. Usually, one varies the choice of scale about the physical 
scale by a factor of two or so, but as discussed earlier, the FAC scale defined in eq. ( p.l2| ) is 
an attractive alternative choice in that the known ultraviolet logarithms are resummed . 
For both of these variables, the FAC scale is significantly less than the physical scale, for 
example, for -Bmm , O.OGa/s. This has the effect of increasing a^, thereby increasing 

the NLO prediction and in both cases, we see much improved agreement at larger values of 
the observable. At smaller values, and particularly in the region where the data turns over 
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Figure 2: The leading order (dashed) and next-to- leading order (solid) predictions evaluated 
at the physical scale fi = y/s = Mz for (a) l/cxhad ■ dcr/dB^i^ compared to the published 
DELPHI data and (b) the difference between data and NLO theory (normalised to 



NLO). The short-dashed line shows the next-to-leading order prediction using the FAC scale 
(see eq. (|2l^)). 



the agreement is still poor. This, of course, is where the infrared logarithms are large and 
need to be resummed. Furthermore, we also expect non-perturbative hadronisation effects 
or power corrections to influence the perturbative shape of the distribution we have been 
concerned with [^, These contributions (together with resummation of the infrared 
logarithms) have played an important role in extracting useful information from analyses of 
three jet shape variables, and are likely to be important in analysing four jet event shapes. 

A similar comparison of the perturbative predictions for the jet transition rates with the 
DELPHI measurement^ is made in Figs. | and ^ Once again, the LO distribution lies well 
below the data. This time, the NLO prediction lies much closer to the data for most of the 
available kinematic range. The FAC scale rate usually lies above the NLO prediction so that 
the data lies within the range of uncertainty engendered by the renormalisation group. 

For completeness. Figs. and ^ show the DELPHI data and perturbative predictions 
for the D parameter and Tminor repectively. As expected from the analysis of Nagy and 
Trocsanyi the LO prediction for D is too low by a factor of about four, while at the 



physical scale fi = ^/s = Mz the NLO distribution is roughly twice as large but still lies 
a factor of two below the data. However, for the FAC scale (which for the D parameter is 
approximately 0.035^^) the prediction overshoots by 50% or so for D > 0.1 where the fixed 
order prediction is least affected by large infrared logs. 

''The DELPHI data gives the differential jet rate rather than the jet transition variable. Up to a small 
(~ few %) correction from five (and more) jet events falling into a four jet configuration, the two quantities 
coincide at fixed order. 
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Figure 3: The leading order (dashed) and next-to- leading order (solid) predictions evaluated 
at the physical scale fi = ^/s = Mz for (a) l/o"had ■ da/d{Ml/ s) compared to the published 
DELPHI data and (b) the difference between data and NLO theory (normalised to 
NLO). The short-dashed line shows the next-to-leading order prediction using the FAC scale 
(see eq. (|2l^)). 



The importance of resumming these logs is clearly shown in Fig. ^ where the Tminor 
distribution is shown. For Tminor > 0.1 the data again lies between the next-to- leading order 
predictions at the physical and FAC scales (which encompass an uncertainty of about a factor 
of almost three for Tminor ~ 0.2). However, the turn-over at Tminor = 0.1 cannot be modelled 
without resumming the large logs which cause the perturbative prediction to grow rapidly. 
The same is true at small values of the light hemisphere mass and narrow jet broadening 
although there the effects are less pronounced because of the choice of bin sizes. 
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Figure 4: The leading order (dashed) and next-to-leading order (solid) predictions evaluated 
at the physical scale fi = y/s = Mz for (a) l/a^ad ■ da/dA compared to the published 
DELPHI data and (b) the difference between data and NLO theory (normalised to 
NLO). The short-dashed line shows the next-to-leading order prediction using the FAC scale 
(see eq. (|2J^)). 
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Figure 5: The leading order (dashed) and next-to- leading order (solid) predictions evaluated 
at the physical scale /i = y/s = Mz for (a) l/ahad ■ dcr/dy( compared to the published 
DELPHI data and (b) the difference between data and NLO theory (normalised to 
NLO). The short-dashed line shows the next-to-leading order prediction using the FAC scale 
(see eq. (^^)). 
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Figure 6: The leading order (dashed) and next-to- leading order (solid) predictions evaluated 
at the physical scale fi = y/s = Mz for (a) l/cxhad • da/dy^ compared to the published 
DELPHI data and (b) the difference between data and NLO theory (normalised to 
NLO). The short-dashed line shows the next-to-leading order prediction using the FAC scale 
(see eq. (|2J^)). 
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Figure 7: The leading order (dashed) and next-to-leading order (solid) predictions evaluated 
at the physical scale = a/s = Mz for (a) l/ahad ■ dcr/dD compared to the pubhshed 
DELPHI data and (b) the difference between data and NLO theory (normalised to 
NLO). The short-dashed line shows the next-to-leading order prediction using the FAC scale 
(see eq. (^^)). 
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Figure 8: The leading order (dashed) and next-to- leading order (solid) predictions evaluated 
at the physical scale /i = ^/s = Mz for (a) l/cThad ■ d<^/dT^mor compared to the published 



DELPHI data and (b) the difference between data and NLO theory (normalised to 
NLO). The short-dashed line shows the next-to-leading order prediction using the FAC scale 
(see eq. (PTT^)). 
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6 Conclusions 



In this paper we have introduced a new Monte Carlo program for the calculation of 4 jet 
like observables in electron-positron annihilation. This program, EERAD2, is based on the 
known squared matrix elements for 7* — > 4 and 5 partons and numerically implements the 
necessary cancellations between the different final states using a hybrid of the commonly 
used subtraction and slicing schemes. This infrared cancellation scheme is detailed in the 
appendix. 

We have checked that the numerical results obtained with EERAD2 are consistent with 
the two other available four jet programs MENLO PARC and DEBRECEN by recalculating the 
distributions for some of the previously determined four jet event variables (such as the D 
parameter, thrust minor and the jet transition rate for the Durham jet algorithm) as well 
as the four jet rate. Within the (estimated) statistical Monte Carlo errors, there is excellent 
agreement. 

We have also presented the leading and next-to-leading order scale independent coeffi- 
cients for some previously uncalculated observables; the light hemisphere mass, narrow jet 
broadening, Aplanarity and the four jet transition variables with respect to the JADE and 
Geneva jet finding algorithms. For each of these observables, the next-to-leading order cor- 
rections calculated at the physical scale significantly increase the rate compared to leading 
order (see fig. |1|). The renormalisation scale dependence is also rather large. 

Furthermore, for each of these observables, we have made a comparison with the published 
data collected at LEP 1 energies by the DELPHI collaboration. With the exception of the 
?/4 distributions, the data lies well above the NLO prediction, apart from in the infrared 
region where the NLO prediction grows rapidly (and where resummation of large logarithms 
is essential). Using the FAC scale, which is of the order of 0.1-1% of the physical scale, 
increases the predicted rate and in general produces a slightly better agreement with the 
data. Taking together the size of the NLO corrections, the renormalisation scale dependence 
and poor agreement with data it appears that the inclusion of even higher order corrections 
and/or power corrections will be necessary to extract any useful information from these 
observations. 

Although at first sight this may seem discouraging, it is instructive to compare this 
situation with the well-established results for the 3 jet-like variable 1 — Thrust. In this 
case, the next-to-leading order coefficient is also large compared to the leading order term, 
resulting in a K factor at the physical scale which varies from 1.4-1.6 throughout most of 
the kinematic range of the variable. These values are very similar to the case for the four 
jet observables that we have already discussed. In addition, it is also possible to compare 
the pure perturbative prediction for 1 — Thrust with the DELPHI data, as we have done for 
the four jet variables in section ^ This yields results which are qualitatively very similar 
to those shown for the D-parameter and Tminor in Figs. |^ and ||. In these cases, it is clear 
that the perturbative prediction - either at the physical scale or using the FAC choice - is 
insufficient to describe the data. However, after resummation of large logarithms and the 
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inclusion of non-perturbative power corrections proportional to the inverse of the energy 
scale of the process, 1/Q, it is known that the data for 1 — Thrust can be described very well 



|37|. In fact the thrust distribution forms a text-book example of how to interpret hadronic 



final states within QCD. 

From this we conclude that discarding four jet-like event shape variables as unreliable 
at next-to-leading order is premature without proper consideration of the types of non- 
perturbative terms that have been successfully included for a variety of three jet-like observ- 



ables. In this we disagree with the conclusions of reference [|TT 



In addition to these correction terms, all of the distributions that we have considered ex- 
hibit the infrared logarithmic behaviour (as the observable tends to zero) that is inevitable 
in a fixed order calculation. These logarithms are present in the coefficients Bq^ and C04 
and can be parametrized in the form shown in eq. ( p.l3 ). By performing a fit to the dis- 



tributions at low enough values of the variable O4 one should be able to extract the values 
of the coefficents Anm and, where possible, perform exponentiation to resum the infrared 
logarithms. 
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A Cancellation of infrared singularities 



In order to compute next-to-leading order quantities in perturbation theory, it is necessary 
to combine the contribution from ?7,-parton one-loop Feynman diagrams with the {n + 1)- 
parton bremstrahlung process. The virtual matrix elements are divergent and contain both 
infrared and ultraviolet singularities. The ultraviolet poles are removed by renormalisation, 
however the soft and collinear infrared poles are only cancelled when the virtual graphs are 
combined with the bremstrahlung process. Although the cancellation of infrared poles can 
be done analytically for simple processes, for complicated processes, it is necessary to resort 
to numerical techniques. 



A.l A simple example 

The numerical problem has been nicely formulated by Kunszt and Soper |^| by means of a 
simple example integral, 

X=lim|^'^x^F(x)-^F(0)|, (A.l) 

where F{x) is a known but complicated function representing the (n + l)-parton brem- 
strahlung matrix elements. Here x represents the angle between two partons or the energy 
of a gluon and the integral over x represents the additional phase space of the extra parton. 
As a; — > 0, the integrand is regularised by the x"^ factor as in dimensional regularisation, 
however, the first term is still divergent as e ^ 0. This divergence is cancelled by the sec- 
ond term - the n-parton one-loop contribution - so that the integral is finite. A variety of 
methods to compute I have been developed. 

The method used by Ellis, Ross and Terrano [|l| is also known as the subtraction method. 
Here, a divergent term is subtracted from the first term and added to the second, 

I = limj/'— x^(F(x)-F(0)) + F(0) x^--F(0)l 
I^Jo X Jo X e j 

- {F{x) - F(0)) , (A.2) 

X 

so that the integral is manifestly finite. This method has the advantages of requiring (in 
principle) no extra theoretical cutoffs and making no approximations. However, in practice, 
there are still large cancellations in the numerator and there is a hidden parameter which 
cuts the integral off at the lower end. Recently, this technique has been developed to describe 



general processes [16|. 



An alternate approach known as the phase space slicing method has also been widely 
used 0. The integration region is divided into two parts, < x < 6 and 6 < x < 1. In the 
first region, the function F{x) can be approximated by F{0) provided the arbitrary cutoff 
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5< 1, 



/■I fj'Tf 

~ / — F(a;) +F(0)ln(5). (A.3) 

Js X 



This method is extremely portable ||T8| since the soft and collinear approximations of the 
matrix elements and phase space are universal. This makes it easy to apply to a wide variety 
of physically interesting processes. However, the disadvantage is the presence of the arbitrary 
cutoff 6. The integral should not depend on 6, and the 6 dependence of the two terms in 
eq. ( [A.3| ) should cancel. Since the approximations are reliable only when 6 is small, this can 



give rise to numerical problems. 

Finally, we can imagine combinations of these two approaches - the hybrid approach. 
There are two scales in the problem, S and A. In the region < x < 6, we adopt the slicing 
procedure, while in the range S < x < A we add and subtract an analytically integrable set 
of universal terms, E{x), to eq. (|A.3|) , 



Xr^ f ^F{x) + F{0)ln{6)- f —E{x)+f —E{x), (A.4) 
Js X Js X Js X 

which on rearrangement yields, 

X ^ /' ^ {F(x) - E(x)e(A - x)) + r —E(x) + F(0) \n(6). (A.5) 
Js X Js X 

Because we explicitly add and subtract the same quantity, there can be no dependence on 
A which can therefore be taken to be large. However, the slicing approximation requires 
5 — s> 0. For this approach to be useful, two conditions must be satisfied. First, the second 
term in eq. ([A.5|) must be evaluated analytically without making any approximation in the 
phase space and should produce a term —F{0) ln(5) from the lower boundary that explicitly 
cancels the third (slicing) term. This allows the limit 5 — to be taken (inasmuch as that 
can be achieved numerically). Second, F{x) ~ E{x) as a; — and more usefully E{x) is 
smooth and as close to F{x) as possible over the whole range of x < A, so that the first term 
in eq. ( [A.5| ) can be safely evaluated numerically. This is the technique we have adopted in 
this paper. 



A. 2 Singular behaviour of Matrix Elements 

Clearly the choice of the subtraction function E{x) requires some care, as does the integration 
over the phase space variables x. To help us do this in a sensible way, we first recall the well 
known singular behaviour of the matrix elements. This is most clearly seen by decomposing 
the amplitude according to the various colour structures. For example, in the process, 

e~^e~ ^ qq + n g, 
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the amplitude can be written as, 

A4(gi, Q2; 1, . . . , n) = S;^\Qi; 1, . . . , n; Q,)Vr 
where the hadronic current is given by, 

P(l,...,n) 



(A.6) 



(A.7) 



Here, 1, . . . , (^2) represents the colourless subamplitude where the gluons (with 

colours ai, . . . , a^) are emitted in an ordered way from the quarks (with colours ci and C2). 
The colour matrices are normalised such that, 

Tr(T'^'T"-') = 

By summing over all permutations of gluon emission, all Feynman diagrams and colour 
structures are accounted for. On squaring, we find that for n > 1 the leading colour piece is 
simply. 



N 



E 

P(l,...,n) 



S^{Q,-X....n-Q^)V^- + 0[^]]. (A.8) 



The subleading terms are slightly more complicated, but may be straightforwardly obtained 
by considering diagrams with gluons replaced by one or more photons. | 

The advantage of the colour decomposition is that the ordered subamplitudes have par- 
ticularly simple singular limits. For example, in the limit where gluon u is soft, we have 
the QED-like factorisation into an eikonal factor multiplied by the colour ordered amplitude 
with gluon u removed, but the ordering of the hard gluons preserved. 



S^{Qi] l,...,a,u,b,...,n; Q2)V^ 

Saubi^Sabj Saui ^ub) '^^(Ql, 1, . . . , 0-, 

with the eikonal factor given by. 



(A.9) 



n-Q,)V^ 



Saubi^abj ^auy ^ub) 



^au^ub 



(A.IO) 



Similarly, in the limit where two particles become coUinear, the sub-amplitudes factorise. For 
example, if gluons a and b become coUinear and form gluon c, then only colour connected 
gluons give a singular contribution. For example. 



, 2 



5^(Qi;l,...,a,6,...,n;Q2)V^^ ^ P„{z, Sab) S^iQi;l, . . . ,c, . . . ,n;Q^)V^ . {A.ll] 



^Here, we have focussed on the two quark process, however, the same type of colour decomposition can 
be appUed to the four quark process, e+e" qqQQ+ {n—2)g, (see for example Ref. |o|). The structure 
appears to be more complicated, however the singular behaviour of individual contributions follows the same 
pattern. 
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For gluons that are not colour connected, there is no singular contribution as Safe — > 0, and, 
when integrated over the small region of phase space where the collinear approximation 
is valid, give a negligible contribution to the cross section. Here z is the fraction of the 
momentum carried by one of the gluons and, after integrating over the azimuthal angle of 
the plane containing the collinear particles with respect to the rest of the hard process, the 
collinear splitting function Pgg^g is given by, 

Ps9^9i^,S) = -^Pgg^giz) (A.12) 

where the usual Altarelli-Parisi splitting kernel with the colour factor removed is given by 
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Similar splitting kernels exist for other combinations of collinear partons 

P,,-.giz) = (z' + il-zf), (A.15) 

with Pgg^q{z) = Pgg^q{l — z). As bcforc, the colour factors have been removed and azimuthal 
averaging of the collinear particle plane is understood. 

In both the soft and collinear limits, the colour ordered squared amplitudes factorise 
into a squared amplitude containing one less parton multiplied by a factor that depends on 
the the unresolved particle and the two adjacent 'hard' particles. We view the two 'hard' 
particles as an antenna from which the unresolved parton is radiated. It therefore makes 
sense to divide the phase space in a similar way and to treat the subtraction term as the 
singular factor for the whole antenna integrated over the unresolved phase space. 



A. 3 Phase space factorisation 

Let us consider an [n + 1) particle phase space described by momenta pi with pf = for 
i = 1, . . . ,n . If the total centre of mass energy is Q, then let us denote the phase space by, 
dPS{Q'^;pi, . . . ,Pn)- As discussed above, we wish to relate the full (n + 1) particle phase 
space to an n particle phase space whenever one of the original (n + 1) particles is unresolved. 
Let the unresolved particle be labelled by u and the two adjacent hard particles by a and b, 
then the phase space can be factorised as, 

dPS{Q'^;pi, . . . ,p„) = dPS{Q'^;pi, . . .,Paub, ■ ■ ■ ,Pn) ^tt^ dPS{saub; Pa,Pu,Pb), (A.16) 

where Paub = Pa + Pu + Pb and = Saub- To factorise the phase space into an n particle 
phase space multiplied by a factor containing integrals over the unresolved invariants Sau 



26 



and Sub that appear in the singular hmits of the matrix elements, we multiply the r.h.s. of 
eq. (Oi by, 

dPS{sAB;PA,PB) I J dPS{sAB;PA,PB), (A.17) 

where particles A and B have momenta pa and pb such that, paub = Pab = Pa + Pb, 
p'a = P% = and Saub = sab- In other words. 



ds 



dPS{Q'';pi, ...,Pn) = dPS{Q^;pi, . . . ,pab, ■ ■ ■,Pn) ^^rr^ dPS{sAB;PA,PB) x dPS' 



smg 



27r 

dPSiQ^;Pu ...,pA,PB,...,Pn)y< dPS^'^^. 



(A. 18) 



As desired, we have the phase space for an final state containing n lightlike particles multi- 
plied by dPS'**™s. Working in four- dimensions and after integration over the Euler angles. 



dPS{Saub;Pa,Pu,Pb) 
JdPS{sAB;PA,PB) 

^Q^2 ^ o-ub^-^ au^-^ ub^-^ ab^ -^au -^ub -^ab) 



(A.19) 



where Xij = Sij/saub- For this to work, a mapping must exist that determines pa and pb for 
a given set of momenta pa, Pb and pu- Many choices are possible ||T7| , |42| and we choose the 
symmetric mapping of f^ . 



PA 

Pb 

where, 
and. 



1 + P 
1-p- 



Sab + 

+ p - 2ri) 



Pa + ripu + 



Sab + 



Pa+{1- ri)pu + 



1-P + 

1 



Sa«(l - p - 2ri 



Pb, 



1 + P 



Sab ~l~ Sub 



Sab ~l~ Sub 



ri 



Sub 



+ Sub 



P 



'^Ib + i^au + Sub)Sab + 4ri(l - ri)SauSub 



Pb, 

(A.20) 
(A.21) 

(A.22) 



SabSaub 

Note that this transformation approaches the singular limits smoothly. For example, as 
Sau 0, then ri 1, p ^ 1 and Pa~^ Pa+ Pu, Pb Pb- 



A. 4 Antenna factorisation of the Matrix elements 

Having f actor ised the phase space, we now wish to find the analogues of the subtraction 
functions E{x) discussed in Appendix |A.1| . These functions should ideally be valid over the 



whole of the antenna phase space dPS^^^^ and, in the soft and collinear regions must match 
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onto the singular limits discussed in Appendix [A.2| . In other words, for a given (n + 1) 
particle amplitude, in the limit where particle u is unresolved. 



S^{. ..,a,u,b,.. .)¥'■ 



A 



aub 



S^i...,A,B,...)V^ 



(A.23) 



where we have replaced the antenna comprising a, u, b by the hard partons A and B to 
obtain an n particle amplitude. The antenna function Aaub depends on the momenta of the 
radiated particles a, b and u, but the n particle amplitude jiS^jV^p does not. 

The leading colour contribution to an observable cross section from an (n + 1) particle 
final state with a particular colour ordering is proportional to. 



71+1 



S^{...,a,u,b,...)V'' J^n+i) dPS{Q'';...,Pa,Pu,Pb,...), (A.24) 



where the observable function j7(„+i) represents the cuts applied to the (n + 1) particle phase 
space to define the observable. Using the factorisation of the matrix elements defined in 
eq. (|A.23|) , when particle u is unresolved we should subtract, 



n+l 



A, 



aub 



S^i...,A,B,...)V'' dPSiQ^;...,pa,Pu,Pb,...] 



(A.25) 



from the (n + 1) particle contribution and, using the phase space factorisation of eq. (|A.18|) , 
add, 



n+l 



Aaub dP5^'"*^ S,{. ..,A,B,.. .)V^ :/•(„) dPS{Q'; . . .,pa,Pb, 



(A.26) 



to the n particle contribution where both the observable function J' and matrix elements 
depend only on the momenta of the n remaining hard partons. Note that for any 
infrared safe observable, in the limit that one particle is unresolved, J^(n+i) J'{n)- In the 
subtraction term eq. ([A.25| ), we use the transformations of eq. ( |A.2CI|) to map the momenta 
Pa, Pu and pb defined in the {n + 1) particle phase space onto the momenta pa and ps used 
in the ra-particle matrix elements and observable functions. In eq. ( |A.26 ), all dependence on 
the momenta of particles a, b and u may be integrated out to give the antenna factor, JF, 



^ab{sab) 



Aauh dPS^ 



(A.27) 



multiplying the n particle cross section (for a given colour ordered amphtude), 
(^1 \S,{...,A,B,...)V^\ J^^)dPS{Q'-...,pA,PB,...). 



(A.28) 



The full set of subtraction terms is obtained by summing over all possible antennae. 

The Dalitz plot for the [AB) — > aub phase space is shown in Fig. |^. In the hybrid scheme 
we are implementing, we use the slicing method of |jl8| in the region mm{sau, Sub) < and 
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Figure 9: The phase space for the decay (AB) auh. The cut min(saM, Sub) = S with 
6 = 0.1 Saub is shown as a sohd hne while min{sau, Sub) = A is shown as a dashed hne for 
A = 0.25 Saub- The region mm{sau, Sub) < S defines where the shcing approach is utihsed, 
with the soft and coUinear regions demarked by dotted hues. Antenna subtraction is apphed 
when 6 < mm{sau, Sub) < A. 



the subtraction scheme in the region, 6 < mm{sau, Sub) < A. In the shcing region, the phase 
space and soft and colhnear approximations to the matrix elements are kept in D = 4 — 2e 
dimensions to regularise the singularities present when either invariant vanishes. Using the 



approach of [0, there are three separate contributions (a) soft gluon when max{sau, Sub) < 5, 
(b) a and u colhnear when Sau < S but Sub > ^ and (c) b and u collinear when Sub < ^ but 

Sau > S. 

Before turning to the explicit forms for the antenna subtraction terms, we note that while 
quarks are only directly colour connected to one particle - a gluon or antiquark, gluons are 
directly connected to two particles - the gluon (or quark) on either side. Therefore,while 
the quark (or antiquark) appear in a single antenna, gluons appear in two. This gives an 
ambiguity in how to assign the collinear singularities of a pair of gluons to each antenna. 
Later we will exploit this ambiguity to make the antenna functions Aaub for different pairs 
of hard partons finite simpler. 



A. 4.1 Quark- Antiquark antenna 

Let us first consider a system containing a quark, antiquark and a gluon. This is produced 
by an antenna comprising of a hard quark and antiquark pair that decays by radiating a 
gluon. Any function that has the correct soft gluon and collinear quark/gluon singularities 
in the appropriate limit is satisfactory. Here the hard particles in the antenna are Q and Q 
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which radiate to form g, q and the gluon g. A suitable choice for the antenna function is, 

\s,{Q-Q)y'? 

( _j_ _j_ '^•^ab-^aub \ 29) 

^aub \ -^ub -^au -^au-^ub / 

Because this is proportional to the three parton matrix elements, \S^{q; g ^\'^, it auto- 
matically contains the correct soft and coUinear limits. Furthermore, it is smooth over the 
whole three particle phase space and singularities only appear in the Sau — ^ and Sub — ^ 
limits. 

Explicitly integrating over the antenna phase space for 6 < min(sa«, Sut) < A we find. 

Since we intend to take the 5 — > limit, the terms of 0{6) may be safely neglected. The S 
independent function is given by, 

^f^\f_ (^) + ^ _ 2Li2(x) +(l-2x + ^]\n( V (A.31) 



27r y V ' ' 2 ^ ' \2 2 \ X 



A. 4. 2 Quark-GIuon antenna 

For antenna made of a quark Q and gluon G, there are two possible ways of radiating. Either 
a gluon can be radiated so that a quark-gluon-gluon system is formed, or the gluon may split 
into a antiquark-quark pair. This latter possibility is subleading in the number of colours 
and the discussion of situations like this is deferred to sec. IA.4.4. 



For a quark-gluon-gluon system there is a less obvious choice of antenna function, par- 
ticularly since the singularity that is produced when the gluon splits sits in more than one 
antenna. If, in the collinear limit, the gluon splits into an unresolved gluon u which carries 
momentum fraction z and a hard gluon b with momentum fraction 1 — z, the antenna func- 
tion should naively be proportional to Pgg^g which is singular as z — > and z 1. This 
corresponds to singularities as both Sub and Sab 0. However, because the collinear 
singularity sits in more than one antenna - the two gluons also occur in a second antenna 
where the role of the two gluons is reversed - we can make use of the A^ = 1 supersymmetric 
identity to rewrite Pgg^g as, 

^gg-^g ~ Pqg-*q + ^gg^i ~ ^qg^g- (A. 32) 

The soft singularities as ^ — > are contained in Pgq~^q while those as 2; ^ 1 are in Pgg^q- We 
therefore divide Pgg-,g amongst the two antennae such that Pgq~,q sits in the antenna where 
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gluon u is unresolved. The z — ^ 1 singularities are placed in the antenna where the role of 
the two gluons is reversed. The remaining Pqq^g may be divided between the two antennae 
according to choice. With a slight modification due to the Pqq^g term, the antenna function 
used for the QQ antenna has the correct limits, so that, 

^ggg ~ "^ggQ I I • (A. 33) 

This is again smooth over the whole three particle phase space with singularities only appear- 
ing in the Sau and Sub limits. In particular, as z 0, the collinear limit matches 
onto the soft limit which would not have been the case if we had divided the soft/collinear 
singularities equally between the two antenna. 

After integrating over the antenna phase space for 6 < min(sau, Sub) < A we find. 



with the 6 independent function J^qq is given by. 



^6 2 6 J \ X ^ 

Antennae containing a gluon and an antiquark are described by, 

Agq = Aggi^^'-^b)^ (A- 36) 

and, 

•^gq(^gq) = ^QciscQ)- (A.37) 

A. 4. 3 GIuon-GIuon antenna 

For antenna comprising only gluons, we repeat this SUSY inspired trick for each of the 
resolved gluons so that. 



A.. = Ag, - — + ■ (A.38) 



Saub \ •^ub'^aub ■^au-^aub , 



Note that Kosower has proposed an antenna factorisation for gluonic processes, 

^Kosower _ ^ j {^aub{Xaub ~ ^gh) + X^fc) \ 
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which, in the u/h colhnear hmit regenerates the full Pgg^g splitting function, as well as the 
soft limits. 



Integration of the antenna function Aggg over the whole of the subtraction region yields. 




A. 4. 4 Antenna where a quark-antiquark pair merge 

There are also configurations when two (or more) colour lines are present, one ending in an 
antiquark the other starting with a quark of the same flavour. Here the matrix elements 
have the form, 



(A.42) 



In the collinear limit, the quark-antiquark pinch the two colour lines together to form a 
single colour line. 



S^{. . .,a,q\q,b, . . .)V^ 



(A.43) 



with Pqq^Gi^jS) given by eqs. ( [A.12| ) and ([A.15|) . There is no soft singularity, nor is there 



any dependence on the type of adjacent parton, a or b. Clearly, the quark-antiquark pair 
can sit in two antennae, (a, q, q) and (g, q, b) and we have some freedom of how to assign the 
singularities to the antennae. There are two obvious choices. Either we divide the singular 
contribution equally over the two antennae, or, we place the part of 



z] m one 



antenna and the {1 — zY part in the other (as we did with the three gluon antenna before). 
While there appears to be no preference, we follow this latter route so that the antenna 
function vanishes as the unresolved particle becomes soft. 



A, 



aqq 



and, 



X 



aq 



aqq \-'^qq-^aqq , 



Aqqh Aaqqiji^aq ^ ■^qbi ■^aqq 



qqbj 



(A.44) 
(A.45) 
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Following this procedure and integrating over the whole of the subtraction region yields, 



and, 

^GaSG,)=J'^G'{sGb). (A.47) 

The factor of Np arises because each of the Np quark flavours may contribute. The 5 
independent function is, 

asNF\ f 2x (\ x^\, (\ — x 



+ 7r-^- 7:-7r ln — • (A.48) 



27r/\ 3 6 9 V6 6/ \ x 



A. 5 Leading colour contribution to e+e 4 jets. 

As a pedagogical example, we consider the leading colour contribution relevant for e^e~ 
4 jets. The sub-leading pieces are similarly calculated but the resulting expressions are 



somewhat lengthy due to the many antennae that are involved |^0[. At leading order in the 
number of colours, only the two quark and n gluon process contributes^, so, at lowest order, 
the cross section is given by, 

da^ _ (27r)5 f - l \ /«siV\' 

X J2 \s^{Qi;Gi,G2;Q2)V^\^ J-(4) dPS{Q^;QiGi,G2,Q2)l2, (A.49) 

P{Gi,G2) 

where X2 is the identical particle factor for the two gluon final state. In practice, the 2! 
permutations precisely cancels the identical particle factor of 1/2!, and it is more convenient 
to keep one particular ordering so that. 



ao s \ N"^ \ 27V 



2 



2 



S^{Qi;Gi,G2;Q2)V^' J^^) dPS(g'; QiGi, G2, Qa)- 

(A.50) 

Similarly, the leading colour contribution from the five parton bremstrahlung process is, 
das (27r)^ f - l\ /a,7V^ 



(To s \ N"^ J \ 2tt 

^The four quark process gives a contribution that is suppressed by a factor of Np /N relative to the leading 
colour contribution. Numerically this is an important contribution and, together with the other subleading 
terms, is included in the numerical results presented earlier. 
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3 



(27r)^ (N^ - l\ (asN 



2tx 



<Sf,{qi]9i,92,93;(l2)V^ J(5) dPS{Q ; qi, 91,92,93,92) ■ 

(A.51^ 



Note that jTs projects the five parton momenta onto the four jet hke observable. Once 
again, we can cancel the identical particle factor I3 = 1/3! against the 3! permutations 
of the gluons, and retain the single permutation indicated. For this colour ordering, three 
antennae will contribute, {qi,9i, 92), {91, 92, 93) and {g2, 93, ^2) where in each case the parton 
in the middle is unresolved. In the first antenna, {'Pqi,'Pgi,'Pg2} ~^ {PQdPGi} according to 
eq. (|A.20|) , the slicing cuts are min(sqigi, Sg-^g^) < 5 and the subtraction occurs over the range 
5 < min(sgigi, s^^gj < Similar transformations and cuts act over the other two antenna. 



A. 5.1 Slicing contribution 



For the five parton matrix elements of eq. ( |A.51| ), the sum of infrared singularities from the 
three antennae in the slicing approach gives a contribution to the four particle final state 



which can be read directly from eq. (3.79) of ref. |T8 



daf^ = /2(Qi;G'i,G'2;Q2)da4^°- 



(A.52) 



Retaining only the leading colour contribution (i.e. dropping the contributions from the 
four quark process proportional to the number of quark fiavours). 



R{Qi; Gi, G2, Q2) 



277 / r(l - e) 
3 /47r/i2\^ 



,2\ f 



+ ^ 



2e\ 6 

271 ~ 



197 

H TT 

18 



2 



6 



e ri 



S 



+ 0(e) + 0(5). 



with (at leading order in the number of colours) 60 = llA^/6 and where the sum runs over 
the pairs of adjacent (colour connected) hard partons, ij = QiGi, G1G2 and G2Q2- 



A. 5. 2 Subtraction term 

Since there are three antennae, we subtract three antennae factors, such that the total 
subtraction term is, 



(27r)^ /A^2 _ i\ .^^jY 



iV2 



27r 



dPS{Q^;qi,gi, 02,93,^2) 
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+ A, 



919293 



+ A, 



9292^2 



'(4) 



(A.53) 



Here, we have used the mappings {Pgi , Pgi , Pc,2 } ~^ {PQdPgi} according to eq. ( |A.20p for the 
first antenna. We recall that the subtraction occurs over the range 5 < min(sgigj, s^^^J < A 
and that the observable function is applied to the momenta for Qi, Gi, g^^ and g2- Similar 
procedures are applied to the other antennae. 

However, we must add these terms back to the four parton contribution. Here it is 
simplest to re-identify each of the four particle momenta with the momenta relevant for 
tree level. In other words, for the first antenna, {Pgi , Pgi , Pga } ~^ {PQdPgi} as before and 
V9-i VGi , Pq2 ~* Pq^ ■ This is safe to do since we integrate over the whole four particle phase 
space. Altogether, we have. 



^ QiGi + ^ G1G2 + ^ G2Q2) '^^4°- 



(A.54) 



A. 5. 3 Virtual contribution 



From ref. J^, the matrix elements for the leading colour one- loop contribution to the qqgg 
final state for this colour ordering can be written in terms of a pole structure in e multiplying 
the lowest order matrix elements and a function Ca that is finite as e — 0, 



Ca{Gi, G2) = V(Qi; Gi, Ga; Q^) S^{Qi- Gi, G2; Q^W'' + Ca{G,, G2). 
The divergent factor V is given by, 

'a,N\ ( V{SQ,G,) ^(^G.Q,) VisG,G2) 



where we have introduced the notation, 

'4VYr2(l-e)r(l + e) 



2 e 



Vis) 



\ -s j r(l-2e) 
In terms of cross sections, we have, 

daj = V{Qi- Gi, G2; Q2)da4^° + da^*^'^'*^ 

with, 



2ti 



(A.55) 



(A.56) 



(A.57) 



(A.58) 



tA{G^,G2) J-(4) dPS{Q^;pQ,,pG,.pG2.PQ;) (A.59) 
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Here La is a finite function given in terms of that logaritlims and dilogarithms tliat arise 
in evaluating the one-loop contributions.!] We have ensured that the kinematic singularity 

structure of La matches that of the tree-level |5^(<5i; Gi, 6*2; ^2)^^ • It can be written 
symbolically as, 

LA = Y.Pi^^)'^i^ (A-60) 

i 

where the coefficients -Pi(s) are rational polynomials of invariants. The finite functions Lj 
are the linear combinations of scalar integrals defined in |]^ which are well-behaved in all 
kinematic limits, so that La is numerically stable. Unfortunately, the analytic expression 
for La is rather lengthy so we do not reproduce it here. 



A. 5. 4 Next-to-leading order cross section 

Assembling the various pieces, and applying coupling constant renormalisation. 



,27r/ V 27r ; V " V 27r y er(l 
the NLO four parton contribution is, 

dar° = ddj + d<^ + daf 

= ir(Qi;Gi,G2;Q2)dor4^° + dar'''"*', (A.62) 

where &oY^ and daj'^'^^*'^ are given by eqs. ( |A.49| ) and (|A.59|) respectively with the replace- 
ment as as(yu)- The factor K is the sum of the divergent one loop factor (eq. ( [A.56| )), the 
slicing factor (eq. ( [A.53|) ) and the subtraction term (eq. (|A.54| )), 



+^ QiGi(sQiGi) + ^ GiGiisGiGi) + ^ G^Q^^^G-i^^ 



a 



27r ; V 18 2 







' A 


"SG1G2/ 




■SG1G2 











Similarly, the five parton leading colour contribution to four jet-like observables is ob- 
tained from eqs. ( |A.51D and (|A.53|) , 



^Note that Ca defined here is a factor of 8 smaller than that in Q. 
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evaluated with the running as{fi)- By construction this is finite as any one particle becomes 
unresolved. In the slicing regions, dcr^^*"* = 0, while the phase space regions over which the 
subtraction terms are applied are implicit in the definition of the antenna functions. 

Note that the four-dimensional limit of all cross sections may be taken with impunity 
now that the singularities have cancelled. Furthermore, there is no dependence in K on the 
slicing parameter 6 which may also be taken as small as desired. The subtraction parameter 
A remains, and both dcr^'"^ and do"^'"^ individually depend on it. However, the sum of both 
contributions is independent of the choice of A. The precise value of A can be made bearing 
in mind the numerical stability and speed of the final computer code. For small A, there 
may be sizeable cancellations between the four and five parton contributions, while for large 
A more CPU time is required to evaluate the subtraction terms. For our numerical results, 
we have taken, 

5 = 10-^ A = 10-\ (A.65) 
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